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Abstract 

In this paper, we consider the dynamical evolution of dark vortex states in the two- 
dimensional defocusing discrete nonlinear Schrodinger model, a model of interest 
both to atomic physics and to nonlinear optics. We find that in a way reminiscent 
of their Id analogs, i.e., of discrete dark solitons, the discrete defocusing vortices 
become unstable past a critical coupling strength and, in the infinite lattice, they 
apparently remain unstable up to the continuum limit where they are restabilized. 
In any finite lattice, stabilization windows of the structures may be observed. Sys- 
tematic tools are offered for the continuation of the states both from the continuum 
and, especially, from the anti-continuum limit. Although the results are mainly 
geared towards the uniform case, we also consider the effect of harmonic trapping 
potentials. 
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1 Introduction 



The study of vortices and their existence, stability and dynamical proper- 
ties has been a central theme of study in the area of Bose-Einstein conden- 
sates (BECs) [TP] . In particular, the remarkable experiments illustrating the 
generation of vortices |3 4j5] and of very robust lattices thereof [6117118] have 
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stirred a tremendous amount of activity in this area in the past few years, 
that has by now been summarized in various reviews and books; see for ex- 
ample Pll0|lll|ll2lll3|ll4j . Much of this activity has been centered around the 
robustness of vortex structures in the context of the mean-field dynamics of the 
BECs (which are controllably accurately described by a nonlinear Schrodinger 
(NLS) equation) in the presence of many of the potentials that are relevant 
to the trapping of atomic BECs including parabolic traps [Tf2] and periodic 
optical lattice ones [ToTTE] . Particularly, the latter context of optical lattice po- 
tentials is quite interesting, as it has been suggested that vortices (for example 
of topological charge 5 = 1) will be unstable when centered at a minimum of 
the lattice potential [17], an instability that it would be interesting to under- 
stand in more detail. 

On the other hand, the BECs in the presence of periodic potentials have 
been argued to be well-approximated by models of the discrete nonlinear 
Schrodinger (DNLS) type (i.e., resembling the finite-difference discretization 
of the continuum equation) [T8|19|20|[2T] . In that regard, to understand the 
existence and stability properties of vortices in the presence of periodic po- 
tentials, it would be interesting to analyze the discrete analog of the relevant 
NLS equation. This is also interesting from a different perspective in this BEC 
context, namely that if finite-difference schemes are employed to analyze the 
properties of the continuum equation, it is useful to be aware of features in- 
troduced by virtue of the discretization. 

However, it should be stressed that this is not a problem of restricted im- 
portance in the context of quantum fluids; it is also of particular interest in 
nonlinear optics where two-dimensional optical waveguide arrays have been 
recently systematically constructed e.g. in fused silica in the form of square 
lattices [22]I25] (and, more recently of even more complex hexagonal lattices 
[21]), whereby discrete solitons can be excited. By analogy to their one- 
dimensional counterparts of discrete dark solitons, which have been created 
in defocusing waveguide arrays with the photovoltaic nonlinearity [25], we ex- 
pect that it should be possible to excite discrete dark vortices in defocusing 
two-dimensional waveguide arrays. An especially interesting feature of dark 
solitons that was observed initially in [25] (see also [27]) is that on-site dis- 
crete dark solitons are stable for sufficiently coarse lattices, but they become 
destabilized beyond a certain coupling strength among adjacent lattice sites 
and remain so until the continuum limit where they are again restabilized (as 
the point spectrum eigenvalue that contributes to the instability becomes zero 
due to the restoration of the translational invariance in the continuum prob- 
lem) [26f27j . It is therefore of interest to examine if the instability mechanisms 
of discrete defocusing vortices are of this same type or are potentially different 
and how the relevant stability picture is modified as a function of the inter-site 
coupling strength. 
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It is this problem of the existence, stability and continuation of the vor- 
tex structures as a function of coupling strength that we examine in the 
present work. We consider, in particular, a two-dimensional discrete non-linear 
Schrodinger equation 

.(hpn^m _ |^ nim |2^ m + e A^ n ,m = 0, (1) 



where A-0 n , m = Vn+i,m + Vv»-i,m + VWn+i + i>n,m-i ~ ¥>n,m is tne discrete 
Laplacian. We study the defocusing case when e > 0. In that case, equation (pP) 
is denoted as discrete Gross-Pitaevskii equation in analogy with its continuum 
counterpart [H2I28] . 

We look for time-periodic solutions with frequency u;. Using the ansatz ip ntm (t) = 
\fu ' <t>n,me~ twi we obtain 

C A0 n , m + (1 - |0„, m | 2 )0„, m = 0, (2) 



where we have set C = e/u. The coupling parameter C > determines 
the strength of discreteness effects. The limit C — > +oo corresponds to the 
continuum (stationary) Gross-Pitaevskii equation: 

S + H^-I^O. ,3, 



The case C — > corresponds to the so-called anti-continuum (AC) limit [29J. 

When equation (j2j) is considered on an infinite lattice Z 2 , we look for solutions 
satisfying |0 n ,m| —* 1 when (n,m) — > oo, for which 0„ im vanishes at one 
lattice site, e.g. at (n,m) = (0,0). Such solutions are denoted as discrete 
vortices, or "dark" vortex solitons. If one trigonometric turn on any path 
Max(|n|, \m\) = p around the vortex center changes the argument of 0„ )m by 
2ttS (S € Z), then the vortex is said to have a topological charge (or vorticity) 
equal to S. 

In this paper we numerically investigate the existence and stability of such 
solutions on a finite lattice of size N x N, N being large; our analysis is per- 
formed as a function of the lattice coupling parameter C and we illustrate how 
to perform relevant continuations both from the continuum, as well as, more 
importantly from the AC limit (section 2). We mainly focus on numerically 
computing vortex solutions with vorticity S = 1 and S = 2 (section 3). In sec- 
tion 4, we also obtain such solutions in the presence of an external harmonic 
trap (the latter is typically present in BEC experiments). Finally, section 5 
presents our conclusions and some future directions of potential interest. 
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2 Numerical method 



We compute vortex solutions of (jSJ) using the Newton method and a contin- 
uation with respect to C. The path-following can be initiated either near the 
continuum limit (for C large) or at the anti-continuum limit C = 0, since 
in both cases one is able to construct a suitable initial guess for the Newton 
method. 



For relatively high C, a suitable initial condition for a vortex with topological 
charge S is obtained with a Pade approximation developed for the continuum 



limit in 



We set 



Pn,r, 



,iSa n , m w h ere 



Pn 



2S+2' 
,m 



V 



n 2 + m z 



(4) 



(ai = 11/32, a 2 = ai/12, b\ = 1/3, see reference [28]). 



ex.. 



arctan(m/n) + for n > 1, 
arctan(m/n) + | for n < — 1, 
| (1 — sign(m)) for n = 0. 



Once a vortex is found for a given C, the solution can be continued by increas- 
ing or decreasing C. Although this method was found to be efficient, it remains 
limited to single vortex solutions having explicit continuum approximations. 
Moreover, when the Newton method is applied to continue these solutions 
near C = 0, the Jacobian matrix becomes ill-conditioned (and non-invertible 
for C = 0) and the iteration does not converge. 



In what follows we introduce a different method having a wider applicability, 
and for which the above mentioned singularity is removed. We consider a finite 
N x N lattice with (n,m) e T = {— M, . . . ,M} 2 (N = 2M + 1), equipped 



with fixed-end boundary conditions given below. We set 



R r 



A6 n 



and note R = (R n , 



m)n,mi 







'n, ml n.m 



. One obtains the equivalent problem 



Rn,m (1 ^n,m 

) + Cf(R,< 



(5) 



Cg(R, 



(6) 



where f(R,9) = Re [ e~ ie A(R e ie ) ] and g(R,9) = Im [e~ ie A(Re ie ) ] can be 
rewritten 
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f(R, 



Rn+l,m COS (9n+l,m @n,m) Rn—l,m COS y-':..m 
~\~Rn,m+l COS (0 n ,ro+l — 9 n , m ) + R nm —\ COS (9 n , 



m. @n— l,m) ^Rn,; 
9n,m— 1 ) ? 



9\Ri")n,m Rn+X,m ^i^l (@n+l,m @n,m) Rn—l,m^>^-ifin,m @n—l,m) 

~~\~Rn,m+l Sill (^n,m+l $n,m) Rn,m—l SU1 (9 n ,rn @n,m— l)- 

Now we divide equation (jH]) by C (this eliminates the above-mentioned degen- 
eracy at C = 0) and consider equation §5§ coupled to 

g(R,e) n:m = 0. (7) 



System (JSJ) , ([Tj) is supplemented by the boundary conditions 



R ntm = 1 for Max(|n|, |m|) = M, (8) 
0n,m = 0™ m for Max(|n|, \m\) = M. (9) 

The prescribed value 9™ m of the angles on the boundary will depend on the 
type of vortex solution we look for. In particular, we use the boundary condi- 
tions 9™ m = Sot n , m for a single vortex with topological charge S centered at 
(n,m) = (0,0). 

For C = 0, a single vortex at (n, m) = (0,0) corresponds to fixing _R 00 = 
and R n ,m = 1 everywhere else. Equation ([7]) yields in that case 

Sin (0 n +l,rn ^n,m) SU1 (9 nrn $n— l,m) 

+ sin (9 ntTn+ i — 9 ntm ) — sin (0„ >m — 9 rijm _ 1 ) = 0, (10) 
(n, m) G T \ { (0, ±1), (±1, 0) , (±M, m) , (n, ±M) } 



supplemented by the four following relations at (n, m) = (0, ±1), (±1,0) 

sin (9 ly ±i - 6> ,±i) - sin (0 o ,±i - 0-i,±i) + sin (0 O)±2 - 0o,±i) = 0, (11) 
sin (0 ±2)O - 0±i iO ) + sin (0 ±M - ±1>o ) - sin (0 ±1)O - ±li _ x ) = 0. (12) 

For a vortex with topological charge S — 1, solutions of f[9l- ([T2"l) are com- 
puted by the Newton method, starting from the initial guess 9 n . m = a n>m . 
The symmetries of the problem allow one to divide by four the size of the 
computational domain. Indeed one can take (n, m) E {0, . . . , M} 2 with the 
boundary conditions 9 0>m = a 0>m , 9 n = a n j0 . Solutions on the whole lattice T 
have the symmetries 

9 n ,-m = vr - 9 n>m [2n], 0_„, m = -9 n<m [2n]. (13) 
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These conditions make (flQl) automatically satisfied at (n,m) = (0,0) (#0,0 
need not being specified). Afterwards, the corresponding solution of (JSJ), ©- 
can be continued to C > by the Newton method, yielding a solution 
0„ im = Rn^e 10 "'" 1 of (see section EJ). For higher topological charges, the 
initial guess 4> n ,m = Rn,m& m can be used to compute a vortex solution of 
(T2J) by the Newton method. This is done in section [3] also for S — 2. All these 
continuations are performed with a 10~ 8 accuracy. 



3 Numerical computation of single vortices 

In this section we analyze the existence and stability of discrete vortices cen- 
tered on a single site, as a function of the coupling strength C for fixed-end 
boundary conditions. The stability of the discrete vortex solitons is studied as- 
suming small perturbations in the form of 5ip m ,n = exp(— it)\p n%m exp(-iXt) + 
g n)TO exp(i\*t)], the onset of instability indicated by the emergence of Im(A) ^ 
0; A in this setting denotes the perturbation eigenfrequency. Note that it is 
sufficient to consider the case u = 1 for stability computations, because this 
case can always be recovered by rescaling time. 

Figure [1] compares the computed angles 9 n>m with respect to the seed angle 
a ntm for fixed-end boundary conditions and iV = 81. The most significant dif- 
ferences arise close to the vortex center. This figure also shows the dependence 
on N of the difference between the angles 9 for a given domain size N and for 
a larger domain of size iV + 10. This is done through \\0^ m — #^ 10 || where 
1 1 ■ 1 1 is the oo-norm, and 9^ m represent the angles at a given lattice size N. 
The main contribution of this norm corresponds to the boundary sites. On 
the other hand, the decrease of this norm as a function of N originates from 
the convergence of the configuration to an asymptotic form. 

Figure [2] shows the complementary norm of the 5 = 1 and S — 2 vortices, 
which is defined as [30J: 

^ = EE(l<M 2 -|0n, m | 2 ) (14) 



with \4>oo\ 2 being the background density; in our case, \(f>oo\ 2 = 1. As it can 
be observed in the figure, vortices with S — 1 and S = 2 can be continued 
for couplings up to 0(1) and presumably for all C0- It should be mentioned 



2 In fact, vortices have been continued at least up to C = 10 without any conver- 
gence problems, and their existence in the continuum limit suggests that it should 
be, in principle, possible to identify such structures for arbitrarily large values of C. 
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Fig. 1. (Left panel) The spatial profile of the difference between the computed angles 
and the seed angles in a 81 x 81 lattice at the AC-limit. (Right panel) Dependence 
°f \\@nm ~ ®nrn°\ loo with respect to the lattice size N . In both cases, the lattice has 
fixed end boundary conditions. 
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Fig. 2. Dependence of the complementary norm on the coupling strength C for 
S = 1 and S = 2. 

in passing that the method has also been successfully used to perform con- 
tinuation in the vicinity of the anti-continuum limit, even for higher charge 
vortices such as S = 3. Notice also that all the considered solutions are "black" 
solitons, i.e., the vortex center has amplitude i?o,o = 0. 

Figures [3] and H] show, for S = 1 and S = 2 vortices, respectively, the profile 
|Vy m | 2 = |0n,m| 2 = Rn mi the angles 9 n ,m, the spectral plane of the stability 
eigenfrequencies and a comparison with the angles a„ jm . In all cases, C = 0.2 
is shown, which corresponds to unstable vortices. 

The vortices with 5 = 1 and S — 2 are, respectively, stable for C < C cr ~ 
0.0395 and C < C cr « 0.0425. This instability, highlighted in the case of the 
S = 1 vortex in Fig. [5] can be rationalized by analogy with the corresponding 
stability calculations in the case of dark solitons [26]. In particular, the relevant 
linearization problem can be written in the form: 
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Fig. 3. Vortex soliton with 5 = 1 and C = 0.2. (Top left panel) density Profile; 
(top right panel) angular dependence; (bottom left panel) spectral plane of stability 
eigenfrequencies [recall that the presence of eigenfrequencies with non-vanishing 
imaginary part denotes instability] ; (bottom right panel) comparison of the vortex 
angles with a n ^ m . 
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(15) 



However, by analogy to the corresponding Id problem, the symmetry and 
the high spatial localization of the localized eigenvector at low coupling ren- 
ders it a good approximation to write for the relevant perturbations that 
^Pn,m ~ — 4j>n,m (and similarly for q), by virtue of which it can be extracted 
that the relevant eigenfrequency is A ~ 1 — AC. This leading order prediction 
(as a function of C) for the internal ( "translational" ) mode frequency is based 
on the anti-symmetry of both the real and the imaginary parts of the vortex 
configuration around its central site, in analogy with the anti-symmetry of 
the on-site dark soliton around its central site in the Id analog of the prob- 
lem [26f27j . This feature (whose continuation to the C — > oo leads to a zero 
frequency mode due to the translational invariance of the underlying contin- 
uum model) is an example of the "negative energy" modes that both dark 
solitons (see e.g., [21] and references therein) and vortices (see e.g. [32]) are 
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Fig. 4. Same as Fig. [3] but for S = 2. 



well-known to possess (due to the fact that, although stationary, they are not 
ground states of the respective Id and 2d systems). 



On the other hand, by analogy to the one dimensional calculation, it is straight- 
forward to compute the dispersion relation characterizing the eigenfrequen- 
cies of the continuous spectrum (using {Pn,m, <?n m) = {-P, Q*} exp[z(A; n n + 
k m m)}, deriving a 2 x 2 homogeneous linear system for P and Q and de- 
manding that its determinant be zero) as extending through the interval 
A G [— V64C 2 + 16C, \/64C 2 + 16(7]. Therefore, the collision of the point spec- 
trum (negative energy) eigenvalue with the band edge of the continuous spec- 
trum yields a prediction for the critical point of C cr ~ (2^/3 — 3)/12 « 0.0387 
in good agreement with the corresponding numerical result above. At C — C cr 
the system experiences a Hamiltonian Hopf bifurcation. In consequence, there 
exists an eigenvalue quartet {A, A*, —A, —A*}. When C increases, a cascade of 
Hopf bifurcations takes place due to the interaction of a localized mode with 
extended modes, as it was observed in one-dimensional dark solitons [26J (see 
also [33], [3l] to illustrate the appearance of this phenomenon in Klein-Gordon 
lattices). This cascade implies the existence of stability windows between in- 
verse Hopf bifurcations and direct Hopf bifurcations. For S = 1 vortices, each 
one of the bifurcations takes place for decreasing |Re(A)| when C grows, and, 
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Fig. 5. Real part of the stability eigenfrequencies for 5 = 1. The panels show 
zooms of two different regions. Dashed lines correspond to the predicted eigenvalues 
A ss 1 - AC and A « V64C 2 + 16C. 
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c c 

Fig. 6. Real part of the stability eigenfrequencies for S = 2. The panels show zooms 
of two different regions. 

in consequence, the bifurcations cease at a given value of C, as |Re(A)| of the 
localized mode is smaller than that of the lowest extended mode frequency 
[however, in the infinite domain limit, this eventual restabilization would not 
take place but for the limit of C — > oo]. This fact is illustrated in Fig. [5J A 
similar plot for the case of the 5 = 2 vortex is shown in Fig. [6l When the 
lattice size tends to infinity (N — > oo), the linear mode band extends from 
zero to infinity and becomes dense; thus, these stabilization windows should 
disappear at this limit. To illustrate this point, we have considered lattices of 
up to 201 x 201 sites for the S — 1 and S = 2 vortices and have shown the 
growth rate of the corresponding instabilities in Fig. [7J The maximum growth 
rate (i.e. the largest imaginary part of the stability eigenfrequencies) takes 
place at C « 0.115 for S = 1 and S = 2 and being Im(A) « 0.0845 (0.0782) 
for S = 1 (S = 2). 
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Fig. 7. Imaginary part of the stability eigenfrequencies for 5 = 1 (left panel) and 
S = 2 (right panel), as a function of the coupling strength C. This corresponds to 
the growth rate of the corresponding instability. 

4 Harmonic Trap 



In this section, we consider the effect of introducing a harmonic trap. Thus, 
Eq. (|2J) is modified: 

C A0„ im + (1 - |0 n>m | 2 - 14,m)0n,m = 0, (16) 

with the parabolic potential of the forrrfl] 

V n , m = -^ 2 < m . (17) 



Fig. [8] shows a typical example of such a discrete vortex structure in the 
presence of an external trapping potential. The method presented in Section [2] 
is again used and converges unhindered by the presence of the magnetic trap. 
Notice that to include the trapping effect of the potential, we only modify 
the initial guess proposed in Section [2] through multiplying it by the so-called 
Thomas- Fermi profile of y / max(0, 1 — V n ,m) [H2] ; the resulting guess converges 
even for small values of C (such as the one used in Fig. [8]). These vortices can 
be continued up to C —>■ oo and will converge to the corresponding continuum 
trapped vortices (for a recent discussion of such vortices in the presence of 
external potentials see e.g. |35j). 

The stability of such structures is also examined in Figs. [9] and HOI The sole 
type of instability observed is an oscillatory one, with alternating windows of 

3 The factor 1/C appears when discretizing the continuum equation given that 
C = 1/h 2 , where h the lattice spacing. In particular, we have used r = \/ x 2 + y 2 = 
hVm 2 + n 2 = hr ntm = r njm /VC. 
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Fig. 8. Vortex soliton with S = 1 and C = 0.5 in a harmonic trap with 0, = 0.1. 
(Left panel) density Profile; (Right panel) angular dependence. 

destabilization and restabilization. However, since the harmonic trap is well- 
known [Hl2] to discretize the spectrum of excitations, these windows of insta- 
bility/restabilization are "true" ones (due to collisions of the "negative energy" 
mode of the vortex with the point spectrum of the background), rather than 
artificial ones (caused by the finite size of the computational domain). In fact, 
in this case, the maximum imaginary part of the eigenvalues does not depend 
on the number of grid points used (provided that the domain "encompasses" 
the harmonically trapped vortex). For high enough C, the charge S = 1 vortex 
is always found to stabilize [32]|35] . It is interesting to also note that although 
the fundamental destabilization scenario indicated by the right panel of Fig. [9] 
has very strong parallels with its untrapped analog, the left panel of the figure 
indicates multiple additional collisions for smaller values of C. The negative 
Krein sign of the translational eigenvalue, discussed previously, suggests that 
these collisions should also result in oscillatory instabilities, although this is 
not discernible in the left column of Fig. [9j A relevant clarification to this 
apparent paradox is provided by Fig. [10] which clearly illustrates that the os- 
cillatory instabilities do indeed arise but, in fact, emerge and disappear (the 
latter through inverse Hopf bifurcations) over very tiny parametric intervals 
of C (and are, thus, apparently invisible over the scale of Fig. [9]). 



5 Conclusions and Future Directions 

In the present paper, we examined the discrete analog of continuum defo- 
cusing vortices which are perhaps the prototypical coherent structure in the 
two-dimensional nonlinear Schrodinger equation. We illustrated how to sys- 
tematically obtain such structures through an appropriate continuation of the 
amplitude and phase profiles from the anti-continuum limit, and also discussed 
how to perform such a continuation from the continuum limit (at least for sin- 
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c c 

Fig. 9. Real part (left panel) and imaginary part (right panel) of the stability eigen- 
frequencies for S = 1 as a function of the coupling strength C and a harmonic trap 
with n = 0.1. 




Fig. 10. Same as Fig.[9]with a zoom around the first bifurcation. From this figure, it 
is clear that there is a Hopf bifurcation that destabilizes the vortex and an inverse 
Hopf. This pair of bifurcations takes place in a small interval of C with a length of 
approximately 3 x 10~ 6 . 

gle core vortices). Such a continuation as a function of the coupling strength 
revealed significant analogies between these defocusing discrete vortices and 
their Id analog of the discrete dark solitons, which are stable from coupling 
C = up to a critical coupling and are subsequently unstable for all higher 
couplings up to C — > oo (when they become restabilized). Something similar 
was observed and quantified in the case of discrete vortices. In addition to the 
most fundamental structures of topological charge S = 1, structures of higher 
charge such as S = 2 were obtained by similar means. 

A natural topic for a more detailed future study arising from the present work 
concerns the understanding of multi-vortex bound states and their stability 
properties, as well as their detailed continuation as a function of the cou- 
pling and eventual disappearance as the coupling becomes sufficiently large. 
Another possible direction would be to examine such defocusing vortices in 
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multi-component models (in analogy e.g., to the bright discrete vortices of 
[36] ; see also references therein). There it would be of interest to study the 
similarities and differences of bound states of the same charge versus ones 
of, say, opposite charges. For these more demanding computations (as well as 
possibly ones associated with the 3d version of the present model [37]), more 
intensive numerical computations will be needed which may be aided by virtue 
of parallel implementation [38] • Such studies are currently in progress and will 
be reported in a future publication. 
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